Dielectric Breakdown in a Mott Insulator: Many-body Schwinger-Landau-Zener 
Mechanism studied with a Generalized Bethe Ansatz 
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The nonadiabatic quantum tunneling picture, which may be called the many-body Schwinger- 
Landau-Zener mechanism, for the dielectric breakdown of Mott insulators in strong electric fields is 
studied in the one-dimensional Hubbard model. The tunneling probability is calculated by a metod 
due to Dykhne-Davis-Pechukas with an analytical continuation of the Bethe-ansatz solution for ex- 
cited states to a non-Hermitian case. A remarkable agreement with the time-dependent density 
matrix renormalization group result is obtained. 
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Among nonequilibrium and nonlinear transport phe- 
nomena in correlated electron systems, dielectric break- 
down (destruction of insulating states due to strong elec- 
tric fields) is one of the most basic. In Mott insulators, 
electrons freeze their motion due to strong repulsive inter- 
action and in equilibrium an introduction of carriers 
in a Mott insulator leads to interesting quantum states 
such as high Tc superconductivity in 2D or Tomonaga- 
Luttinger liquids in ID. Now, it is an intriguing problem 
to ask how nonequilibrium carriers behave when elec- 
trons in a Mott insulator start to move in strong enough 
electric fields. 

The nonequilibrium phase transition from Mott insu- 
lators to metals by electric fields has been studied in the 
condensed-matter physics [3, [3, 0, [H]- More recently, the 
problem is attracting interest in the cold atom physics, 
where novel realization of the Mott insulator has been 
achieved in bosonicjl, 0, [l] as well as in fermionic[l0| 
systems. The many-body Landau-Zener mechanism for 
dielectric breakdown has been proposed for fermionic 
systems in ref. [1], and for bosonic systems in ref. Q. 
The correspondence between the Landau-Zener mecha- 
nism and the Schwinger mechanism [ij in strong- field 
QED as well as the relation between the Heisenberg- 
Euler effective Lagrangian and the nonadiabatic geomet- 
ric phase was given in ref. Q (see also ref. In 
ref. [5| an extensive numerical calculation was performed 
to obtain the electric-field induced nonequilibrium phase 
diagram. One important prediction of the Schwinger- 
Landau-Zener picture is that the threshold electric field 
£^th for the breakdown is related to the charge gap A(C/) 
as C-Eth oc A^(J7), which is much smaller than a naive 
guess of ei?th ~ U, i.e., the energy offset between neigh- 
boring sites in a tilted potential [U: the onsite repul- 
sion) . Such lowering of the threshold was experimentally 
observed by Taguchi et al. 0] who measured the I-E char- 
acteristics in a one-dimensional Mott insulator, where a 
quantum origin of the breakdown was suggested from 
a threshold that remains finite in the zero-temperature 
limit. In cold atoms, the effect of the potential gradient 
was studied to probe the excitation spectrum (they 
use the relation eEth ~ C/ to interpret their results). 

However, the Schwinger-Landau-Zener theories have a 



snag in many-body systems: As explained below eqn. 
the Landau-Zener threshold contains a factor that de- 
pends on the system size and diverges in the thermody- 
namic limit, i.e., no breakdown would take place in bulk 
systems, which contradicts with intuition. The purpose 
of the present paper is to resolve this puzzle, where an 
analytic expression for the threshold field strength valid 
in the thermodynamic limit is presented. This has been 
achieved by deriving the quantum transition probability 
utilizing a method due to Dykhne-Davis-Pechukas (DDP) 
formalism which enables us to treat quantum tunneling 
beyond the Landau-Zener picture [Tsl. 

The present approach has another virtue: Besides the 
quantum tunneling approach, there is a non-Hermitian 
approach studied by Fukui and KawakamiQ, where the 
authors incorporated phenomenologically the effect of 
electric fields as differing left and right hopping terms 
(for non-Hermitian models see also [H, [11]). However, 
the relation to experiments was not too clear, since a di- 
rect connection between the ratio of the left- and right- 
going hoppings with the applied field strength was not 
given. In the present derivation, the non-Hermitian for- 
malism emerges naturally, and the two apparently un- 
related theories (i.e., Schwinger-Landau-Zener and non- 
Hermitian) are shown to be in fact intimately related. 
Indeed, the transition probability in the DDP is calcu- 
lated with an analytic continuation of the solution of the 
time-dependent Hamiltonian onto a complex time, and 
the Hubbard model in an electric field is mapped onto a 
non-Hermitian model. In order to complete the calcula- 
tion, we need the information on excited states. This has 
been achieved here for the ID Hubbard model with a non- 
Hermitian generalization of the Bethe-ansatz 1^ T§| ex- 
cited states, i.e., the string solutions [H, [10, HlOir The 
present result turns out to agree with the time-dependent 
density matrix renormalization group result [^ with a re- 
markable accuracy. 

Here we consider the time evolution of electrons in a 
strong electric field E for the one-dimensional Hubbard 
model. 
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FIG. 1: (color online) Many-body energy levels against the 
complex AB flux $ for a finite, half-filled ID Hubbard model 
{L — 10, — = 5, U = 0.5). Only charge excitations are 
plotted. Quantum tunneling occurs between the groundstate 
(labeled as n = 0) and a low-lying excited state (n = 1) as the 
flux $(t) — Ft increases on the real axis, while the tunneling 
is absent for the states plotted as dashed lines. The wavy 
lines starting from the singular points (x) at "l?(t*) represent 
the branch cuts for different Riemann surfaces, along which 
the solutions n — and n = 1 are connected. In the DDP ap- 
proach, the tunneling factor is calculated from the dynamical 
phase associated with adiabatic time evolution (DDP path) 
that encircles a gap-closing point at $(**) on the complex $ 
plane. 



where the electric field is introduced by a time-dependent 
phase <J>(t) = Ft with F ^ eE switched on at t = 0. 
This is one obvious way of introducing an electric field 
through Faraday's law. We have taken the absolute 
value of the hopping as the unit of energy. We study 
a half- filled, nonmagnetic case with numbers of elec- 
trons A^i = Ni = L/2 with L the total number of 
sites. The Mott-insulator groundstate becomes unsta- 
ble when the electric field becomes strong enough, for 
which charge excitations take place due to nonadiabatic 
quantum tunneling [J. In order to describe the pro- 
cess we introduce the adiabatic levels \ipn{^)) that satisfy 
i/($)|V'n(*)) = -B„($)|V'n(*)) with n = 0, 1, . . ., where 
n = corresponds to the groundstate. We neglect spin 
excitations to concentrate on charge excitations. The 
time evolution for t > is described by the time depen- 
dent Schrodinger equation, i-^\ip{t)) = H{t)\ip{t)), with 
initial state \i;{t = 0)) = \iIjo{<^ = 0)). Figure [D plots the 
adiabatic energy levels obtained by exact diagonalization 
for a small system. Nonadiabatic quantum tunneling 
between the groundstate and the lowest charge-excited 
state is most relevant (while the transition to the state 
represented by dashed lines is absent due to symmetry 
reasons). The adiabatic levels are periodic in $ with a 
period 27r/L, so that the tunneling from the groundstate 
to the excited state repeatedly occurs with a time inter- 
val T = 27t/FL. We define the tunneling factor between 
the two states by 70^1 which is related to the transi- 



tion probability P = e '"'-•i for a single tunneling event. 
The solution of the time-dependent Schrodinger equa- 
tion behaves as \tpimT)) - (1 - e"''"-! )™/2e'"(*) |V'o(0) 
with a phase factor a, and the groundstate decay rate F 



defined by |(^.o(<i>(t))|V'(t)>l' 



-r* becomes0 T/L 



g-70^1)^ A naive estimate for the tunneling 
factor can be made by approximating the Hamiltonian 
in the vicinity of the transition by a Landau-Zencr form, 
j^LZ _ which leads to a threshold behavior 

with threshold F^^^ given by [1, \^ 
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where A is the charge gap (Mott gap) and v is the 
slope of the adiabatic levels (w ~ 2 when U is small 
and the system size is small). However, this expression 
should fail when the system size exceeds the localization 
length [i^], since the slope vanish u — > and the levels 
become flat against $. Then, the transition probability 
also vanishes. But this obviously contradicts with a phys- 
ical intuition that dielectric breakdown should take place 
in infinite systems. The point is that quantum tunneling 
can take place even when the levels are flat [23| . 

In order to resolve this problem we introduce the DDP 
method which accommodates the thermodynamic limit 
as we shall see. In the general formalism of DDP the so- 
lution of the Schrodinger equation is extended to complex 
time; The tunneling process is described by an adiabatic 
evolution of the wave function along a path in the com- 
plex plane (DDP path in Fig. [l] displayed for a finite 
system for clarity). The DDP path encircles the point 
t* (exceptional point) on the complex t plane at which 
the two energy levels cross, i.e., i?i($(i*)) = £'o($(t*)). 
There is a branch cut starting from t* at which the two 
Riemann surfaces corresponding to Eq and Ei merge, 
and along a path encirchng t* the solution |-0o) is de- 
formed into the excited state oc |'0i) with a proportion- 
ality factor determined by the complex dynamical phase. 
This gives a DDP tunneling probability P ~ e^'>'o^i with 

[a, [3 am 



o^DP = 2Im5o,i/ft, 
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where S is the dynamical phase given by 



dt'[E,mt')) - Eom'))i 



(3) 



(4) 



with to the starting point on the real axis. 

We want to apply the DDP method (eqns. Q) to 
the Hubbard model, which means that we have to ana- 
lytically continue the solutions to complex $ for the first 
excited state {Ei) as well as for the groundstate (Eq). 
The Hubbard model with the phase factor (eqn. ([T])) is 
exactly solvable with the Bethe ansatz method (see for 
example [1^). This remains the case, for the ground- 
state, even when $ is complex as demonstrated by Fukui 
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and Kawakami However, we have to extend the pro- 
cedure to the excited states (Fig. [5]), which is feasible with 
Woynarovich's method , where our goal is to calculate 
the energy difference — Eo{^) for complex $ and 

perform the integral along the DDP path. The DDP path 
(Figd]) for the Hubbard model starts from <I>o = tt/L and 
ends at = Tr/L + i^cr, where 'i'cr is the value at which 
the gap closes [1]. In the large L limit the path lies on 
the imaginary axis. 

We start with the Lieb-Wu Bethe-ansatz equation for 
an L-site Hubbard model with an imaginary $ = zVP, 



JV, 



Lkj = 2ttIj + iL^ - ^ 6* (sin kj - A^), (5) 

,(6) 



N, 



/3=1 



where kj (A^) are the charge (spin) rapidities, 0{x) = 
— 2arctan(a::/M) with u = U/{At) is the two-body phase 
shift, and = NJ2 (mod 1), ,/„ = {N - + 
l)/2 (mod 1). 

In the infinite-size limit, the Lieb-Wu equation for a 
finite ^ can be solved with the analytically continued 
charge and spin distribution functions [Sj] . If we introduce 
the counting functions Zc{kj) ~ Ij/L and Zs(Aq,) = Jq/L, 
the Lieb-Wu equation in the bulk limit reads 



^ ^ 27r 27r 27r 



(iA6'(sinfc- A)a*(A), (7) 



^.(A) 



1 

'2^ 



dx'e 



sin fc — A)p* (fc) 
A- A' 



a(A'), 



IS \ 2 

where the distribution functions arc defined by p{k) 



(8) 



(a) ImO=0 



(b) ImO^pO 




FIG. 2: (color online) Schematic configurations (displayed 
here for U = 4.0) of the charge rapidities for the lowest charge- 
excited state |i/)i(i\l/)) for »]/ = 0(a), and for a finite ^'(b). 
C corresponds to the groundstate continuum with occupied 
states (green circles, printed grey). In the excited state, two 
holes ki , km (open circles) near n + ib appear in the contin- 
uum, while two rapidities ki , K2 outside the continuum C are 
occupied. The surface represents the real part of the excita- 
tion energy Re£(fc) (plotted here for Ree(fc) > 0, ImA: > 0) 
which gives the energies of holes at ki, k,n- At = 0, ki, K2 
sit on the Ree(fc) = curve. 



dkZc{k), cr(A) = d\Zs{X), and iT*,a,p* are explained 
around eqn. (jlip below. The contours C and 5, 
i.e., the continuum limit of the charge (C) and spin 
(S) rapidities' positions, are of great importance. 
In fact, for the groundstate, the paths are deter- 
mined such that the conventional solutionflTj. po{k) = 



27T 27r ^^'^ Jo cosh ua; 
Jq{uj) cos a; a 



Jo {(jj) cos {lu sink) duj and o'o(A) 



1 roo ^iH"; Losa;A ^ with J„ Besscl's function, extended 

Zir JO cosn'(iuj ' 

to complex k and A solves eqns. (O, (|H1). This determines 
the end point of contour C, which we denote ±Tr + ib (Fig. 
[2]), where b is an increasing function of satisfying Q 



dXe{X + isinh6)(To(A). 



(9) 



We denote the end point corresponding to = 'I'cr to be 

= fecr- The end point of 5 is A = ±7r. 

Woynarovich's construction [13] of charge excitations 
can be applied to the non-Hermitian case (^f ^ 0) with 
the same contours C, S as in the groundstate. The idea 
is to remove two charge rapidities ki , km from C and one 
spin rapidity Ajv^/2 from S to place them on the complex 
k and A planes at positions ki, K2 and A, respectively 
(Fig. [2]) in such a way that the Lieb-Wu equation is sat- 
isfied, which yields 

sin(Ki^2) = A ± m, A = (sin/c; + sin/s,„)/2. (10) 

With these parameters, the Lieb-Wu equation ([5]) for 
charge excitations can be solved by 



^(A) = MX) 



^ i 

LU 1^ cosh[(A — sin fc;)7r/2it] 
1 



cosh[(A — sin fc„i)7r/2t/] 



pik) = po{k) 



1 



■ cos k- 



(11) 



-{cos[Li;(sin fc — sin/c/)] 
+ cos[u! {sin k — sin km)]} duj, (12) 



2nL ""'^ u2 + (sin k - A)^ 
cosfc f°° e^"" 



2'kL ./g coshwu 



with which we can define cf*{X) ~ o'(A) + {1/L)5{X — 
A), p*{k) ^ p{k) - {l/L)6{k ~ ki) - {l/L)6{k ~ km) 
appearing above. Wc note that these equations are 
identical with Woynarovich's, which is natural since the 
operations dk,d\ do not pick up ^P, while controls 
the integration path via eq.®. The energy of the ex- 
cited state can be calculated from p*(fc), which gives 
Ei{^)-Ea{^) = e{ki) + e{km) with i^o the groundstate 
energy, and the s{k) given as 



e{k) = 2u + 2cos(fc) 



(13) 



g uj cosh uuj 



Ji (uj) cos(ti' sin k)duj. 



The lowest excited state is given by setting ki, km — ir+ib 
in the above solution (Fig[2]). We can specify the defor- 
mation of the Bcthe ansatz solution along the DDP path 
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FIG. 3: (a) The threshold field strength F against U obtained 
by the present DDP formalism (solid line; eqn. (IMf) ). Inset; 
The energy difference between the groundstate and the ex- 
cited state against h for various values of U . (b) The decay 
rate V of the groundstate against the electric field F obtained 
by the DDP formalism (solid line; egn. pSf) ). In (a) and (b), 
the symbols represent the time-dependent DMRG result 



(Fig. [T|) as follows. As ^' becomes finite, the end points 
of C, i.e., ±7r -I- ih, move along the imaginary axis mitil 
h reaches &cr at which the gap closes, i.e., Ei — Eq = Q 
(Fig.O inset) i. Meanwhile, Im ki and Im K2 increase 
with where k,2 in particular touches the real axis at 
the critical point. From the DDP formula (eqns. ([3]), ([3])), 
the quantum tunneling probability P = e"'^^"! has 
the threshold electric field, 



F^^^^^ f\E,-E,) 

71" Jo 



-db 



sinh Li 



1 — cosh h 



— cosh b 



„Lj sinh b 



Jiiio) 



^ c^(l + e2"l-l) 



dcu 



Jo {lo) cosh([j sinh b) 



1 



db. (14) 



Its J7-dependence is plotted e in Fig. [3] (a) (solid line), 
which confirms the collective nature of the breakdown 
(i.e., the threshold much smaller than a naive U). In 
other words, the tunneling takes place not between neigh- 
boring sites, but over an extended region due to a leakage 
of the many-body wave function, where the size is roughly 
the localization length . 

Let us now compare the present analytical result with 
the numerical one in Fig. [3] (a), which plots i^jh'^^ 



along with the threshold obtained by the time-dependent 
density-matrix renormalization group (DMRG) for an 
L = 50, open Hubbard chain The agreement between 
the analytical and numerical results is excellent. 

Finally, let us say a few words about the dynamics that 
takes place after the electric field exceeds the threshold. 
There are infinitely many excited states whose energies 
are larger but near |V'i}'s, and tunneling becomes also ac- 
tivated to these states. The net tunneling to such states is 
incorporated in the groundstate decay rate T/L (defined 
above eqn. ([2])). This quantity has been numerically cal- 
culated with the time-dependent DMRG in ref . Q , where 
the single-tunneling formula reduced by an empirical fac- 
tor a < 1, 



T/L 



aF 
'2^ 



ln[l — exp( 



is found to describe the numerical result. The present 
DDP result again exhibits an excellent agreement with 
the numerical one (Fig. [S^b)). This implies that the 
tunneling to higher excited states do not change the 
threshold, while the decay rate is reduced due to the 
pair-annihilation processesQ- We note that the decay 
rate is an experimental observable which can be obtained 
from the delay time of the current (the production rate 
in ref. Fig. 4), and the present theory is consistent 
with the experimental result. The nature of the nonequi- 
librium steady state above the threshold is an interest- 
ing problem, which will be addressed elsewhere where an 
electron avalanche effect is evoked for the metallization. 

In conclusion, we have shown that the DDP theory of 
quantum tunneling combined with a generalized Bethe 
ansatz describes the nonlinear transport and dielectric 
breakdown of the ID Mott insulator. This is the first an- 
alytical result obtained on nonequilibrium properties in 
correlated electron system, and the DDP method is ex- 
pected to have potential applicability to many other mod- 
els and problems. We wish to thank Mitsuhiro Arikawa, 
Yasuhiro Hatsugai and Takahiro Fukui for fruitful dis- 
cussions, and Seiji Miyashita for bringing our attention 
to [2^. HA was supported by a Grant-in- Aid for Sci- 
entific Research on Priority Area "Anomalous quantum 
materials", TO by a Grant-in- Aid for Yoimg Scientists 
(B) from MEXT. 
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